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Abstract 

A cheapest stable nonconforming finite element method is presented for solv¬ 
ing the incompressible flow in a square cavity without smoothing the corner 
singularities. The stable cheapest nonconforming finite element pair based on 
Pi X Pq on rectangular meshes [28] is employed with a minimal modification of 

the discontinuous Dirichlet data on the top boundary, where is the finite 
element space of piecewise constant pressures with the globally one-dimensional 
checker-board pattern subspace eliminated. The proposed Stokes elements have 
the least number of degrees of freedom compared to those of known stable Stokes 
elements. Three accuracy indications for our elements are analyzed and numer¬ 
ically verihed. Also, various numerous computational results obtained by using 
our proposed element show excellent accuracy. 

Keywords: Nonconforming finite element method; incompressible 
Navier-Stokes equations; lid driven cavity problem 


1. Introduction 

The lid driven square cavity has been one of the most popular benchmark 
problems for new numerical methods for the incompressible Navier-Stokes equa¬ 
tions in terms of accuracy, numerical efficiency and so on. To refer only few see 
[4, 8, 17, 18], for instance, and the references therein. The presence of singular¬ 
ities at the upper corners of the cavity is the source of numerical difficulties for 
solving the cavity flow problem. It is usually erroneous to use high-order meth¬ 
ods without handling the corner singularities due to the Gibbs phenomenon. 
Many studies have been carried out to overcome this difficulty. Barragy and 
Carey [6] used a p-version finite element formulation (p > 6) combined with a 
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strongly graded and refined mesh to handle the corner singularities. Other stud¬ 
ies change the boundary condition to overcome this difficulty: see, for instance, 
[20, 21, 32, 33], and the references therein. The latter approach are coined as 
the so-called regularized lid driven cavity problem. The constant boundary con¬ 
dition for velocity is replaced by a function that vanishes at the upper corners 
of cavity [20, 33]. Botella and Peyret [8] solved a regularized cavity problem by 
using a subtraction method of the leading terms from the asymptotic expansion 
of the solution of the Navier-Stokes equations in the vicinity of the corners, 
where the velocity is discontinuous. Sahin and Owens [32] inserted leaks across 
the heights of the finite volumes at the corners between the lid and the verti¬ 
cal walls to handle the corner singularities. Many studies reported that in the 
critical Reynolds number range [7000,8500] Hopf bifurcations occur for the lid 
driven square cavity problem [4, 18, 20, 33]. Bruneau and Saad [9] revisited 
the issue of bifurcation using third-order time discretization schemes with the 
5000 X 5000 finite difference spatial discretizations. They observed the first bi¬ 
furcation occurs between Re = 8000 and Re = 8050. Guermond and Minev 
[24] reported three-dimensional benchmark solutions using a direction splitting 
method introduced in [22, 23]. They also provided two-dimensional solutions, 
which are correct up to at least three digits, for Re = 1000 using the uni¬ 
form 5000 X 5000 MAC stencil. Instead of the square domain, Glowinski et al. 
[21] considered a semi-circular cavity-driven flow with a special time-dependent 
regularization on the Dirichlet data at the two corners: they observed Hopf 
bifurcations around Re = 6600, which is smaller than the case of square do¬ 
main, using an iso-parametric variant of the Bercovier-Pironneau element [7] 
introduced in [20]. 

The purpose of the current paper is to try to solve the lid driven square 
cavity problem without any regularization at the corners, employing noncon¬ 
forming finite element pairs whose degrees of freedom and implementation are 
as cheap as possible. As the nonconforming elements use the values at the 
midpoints of edges as DOFs, instead of those at the vertices, the discontinu¬ 
ity singularities at the corners are naturally treated without any regularization. 
Our nonconforming finite element pairs are based on the two stable noncon¬ 
forming finite element pairs on uniform square meshes [28] introduced for the 
stationary incompressible Stokes problem. The two pairs are briefly described 
as follows: The first of them uses the Pi-nonconforming quadrilateral element 
[30] for the approximation of the velocity field, componentwise, while the pres¬ 
sure is approximated by a subspace of the piecewise constant functions whose 
dimension is two less than the number of squares in the mesh. The second of 
them is a one-dimensional modification of the above finite element pairs to both 
velocity and pressure spaces: the velocity space is enriched by a globally one¬ 
dimensional DSSY(Douglas-Santos-Sheen-Ye)-type bubble function [11, 15, 26] 
while the pressure space is the subspace of the piecewise constant functions 
whose dimension is one less than the number of squares in the mesh in order 
to fulfill the mean-zero property. The stability and optimal convergence re¬ 
sults for these element pairs applied to the stationary Stokes equations with the 
homogeneous Dirichlet boundary condition can be found in [28]. 
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In order to treat the inhomogeneous lid-driven Dirichlet boundary condition, 
we modified the above elements [28] as follows. The boundary condition on 
the interior of the top boundary is handled as usual, but the corner boundary 
condition is specially treated at the two end elements on the top by adding 
two local DSSY-type bubble functions whose values at the midpoints of top 
boundary parts to be (1,0)* and at the midpoint of the other boundary parts to 
be (0,0)*. Indeed, since nonconforming finite element methods can avoid vertex 
values degrees of freedom, the boundary values at the top left and right corners 
are not required. Thus, one can solve the driven cavity problem without any 
regularization of the boundary condition (3). 

We note that the above modified finite elements have the smallest DOFs 
and are easiest to implement among the finite element space pairs that fetch 
all non-spurious piecewise constant pressure fields. Moreover, our finite element 
methods yield nearly divergence free velocity fields. Indeed, |V ■ v/ijdx = 
which is good indication of numerical solver. The factor arises 

from the inhomogeneous boundary data (the finite element pairs introduced 
in [28] for the homogeneous boundary condition yield exactly divergence free 
velocity approximation.) Another indication of superiority of our element is that 
our methods gives substantially smaller volumetric flow rates across horizontal 
and vertical line sections [5] than other methods by a factor of two. They are 
reported in §5. 

The plan of our presentation of this paper is as follows. In the next section 
the lid-driven cavity problem is briefly described. With a brief review on the Pi- 
nonconforming quadrilateral element, a detailed description and implementation 
of our finite element methods are given in §3. Three accuracy indications of our 
numerical solutions are analyzed in §4. Some numerical results are presented 
§5 with comparison to the results of other methods. The last section concludes 
our presentation. 

2. Problem formulation 

Let n = (0,1)^ be the square cavity. Consider the steady-state incompress¬ 
ible Navier-Stokes equations in dimensionless form: 

—-I- (u ■ V)u -f Vp = 0 in II, 

V • u = 0 in n, 

with the Dirichlet boundary condition 

u = g on r with J iz • g ds = 0. 

Here, u and p denote the flow velocity and pressure, v the fluid kinetic viscosity, 
r the boundary of SI, and u the unit outward normal vector to SI. Here, and in 
what follows, bold faces will denote the two-dimensional vectors, functions, and 


( 1 ) 

( 2 ) 
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function spaces. For the driven cavity problem, suppose that the Dirichlet data 
is given by 

{ (1,0)* if 0 < a; < 1 and y = 1, 
arbitrary if a; = 0,1 and y = 1, (3) 

0 elsewhere on 

Notice that the regularity of the boundary value of velocity field: g € H5“'^(r) 
for arbitrary e > 0, which limits the regularity of the solution u G at 

best. The possible highest regularity of the solutions is 

(u,p) G X W°’^{n)/R for all r G (1,2). 

The Sobolev embedding theorem implies (u,p) ^ H^(fl) xL^(n)/R, but (u,p) G 
H^“'^(f2) X for arbitrary small e > 0. See [12] for more details of 

analysis in the case of driven cavity Stokes equations. 

3. A cheapest nonconforming finite element method 

In this section we will begin with a brief review on the Pi-nonconforming 
quadrilateral element [2, 3, 29, 30] Then the stable cheapest finite element pairs 
[28] for the incompressible Stokes equations with homogeneous boundary condi¬ 
tion will be described. In the third part of this section describes the treatment 
of nonhomogeneous boundary condition for the lid-driven cavity problem. Es¬ 
pecially the corner singularities will be taken care of. 

3.1. The Pi-nonconforming quadrilateral element space 

In this paper, we consider the unit square domain = (0,1)^ with uniform 
square meshes. Let {3^h)Q<h<i be a family of partitions of il into Nq = 
disjoint squares Qjfc ofsize/ix/i, h = 1/A, with barycenter ((/ — ^)h,{k — ^)h) , 
for /, fc = 1, • • • , A. We assume that A is an even integer. By A* denote the 
number of interior vertices Vjk = {jh, kh) in so that A/ = (A — 1)^. Set 

^ g L^(f2) \v\Q.f^ G Pi{Qjk) yQjk & ^h, is continuous at 
the mid point of each interior edge in and v vanishes at 
the mid point of each boundary edge in 5^}, 

The global basis functions of can be defined vertex-wise: for each interior 

vertex Vjk in 5^, define (fjk G such that it has value 1 at the midpoint of 

each interior edge whose end points contains the vertex Vjk and value 0 at the 
midpoint of every other edge in 5^. Then the Pi-nonconforming quadrilateral 
element space [29, 30] is given by 

^ I ^ ^ ^ ajk4>jk I ctjk e R Vj, fc i , dim(.^”^) = A*. 

[ hk=l ] 
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3.2. The DSSY-type finite element space 

The DSSY nonconforming element space on a reference domain Q := [—1, 1]^, 
with vertices Xi = (1,0), X 2 = (0,1),X3 = (—1,0),X4 = (0, —1), is defined by 

DSSY{Q) = Span{l,x,y, 6»fc(x) - Ok{y)}, 


where 

£ = 0 , 

Si{t) = \t^- i = 1, 

[t^-ft^ + lt^, i = 2. 

The reference DSSY basis functions have the form 

= \ + \^- - Oiiy)), 

Y(x) = l + ly + Imx) - OM), 

= i - ij- ^{di{x) - Oiiy)), 

= i - iy + - OM), 


such that (xk) = Sjk, the Kronecker delta. In what follows, we fix £ = 1. 

Let Fq : Q — >■ Q be a bijective affine transformation from the reference 
domain onto a rectangle Q. Then DSSY{Q) is defined by 


DSSY{Q) = {u o 


veDSSYiQ)'^ . 


( 4 ) 


Then the DSSY-type finite element space [10, 11, 15, 26] is defined by 

DSSY^^ = {ve L^ip.) I v\q G DSSY{Q) VQ G 

V is continuous at the midpoint of each interior edge 

and vanishes at the midpoint of each boundary edge in ^}. 


Remark 3.1. For i = 0, the DSSY-type nonconforming element, or the velocity 
components in the CDY(Cai-Douglas-Ye) Stokes element, is identical to the 
rotated Qi element of Rannacher and Turek [31 j. The difference between and the 
DSSY-type nonconforming elements (with £ = 2,3^ and the rotated Qi element 
is that the former satisfies the mean value property u ds = vfme) on each 
edge e in where denotes the midpoint of e. See [26] for more details. 


3.3. The stable cheapest finite element pairs: homogeneous Dirichlet boundary 
case 
Set 


^0 = {Ph^ 


= TjkXQik I 7jfc e R; [ dx = 0 i C Llp), 
J,k=l I 
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where XQjk denotes the usual characteristic function. Denote by the sub¬ 
space of by removing the globally one-dimensional global checkerboard pat¬ 
tern from One way of forming the basis for the {Nq — 2)-dimensional 

space can be described as follows. Let D = U be a decomposition 
of fl into the disjoint unions of red and black rectangles = Uq^^rQ and 
= Uq^^bQ, where 


= {Qjk e I j -I- fc is an even integer}, 

= {Qjk & ■S^h \ j -|- fc is an odd integer}. 

We are now in a position to form the two -dimensional subspaces 

of as follows: 


^R.o 


7>h 

B,0 


= \ph 


= SPh 


= Y] ^jkXQjkDQR I Xjk e K; / dx = 0 > C Ll{n), 

j,k=l ) 

= Y TifeXQ.fcnoB I Ijk eR; J^Ph dx = 0 > C Lo(^)- 


3,k=l 


Then it turns out that q, from which the basis functions 

for is built in a standard way by taking the union of the basis functions of 
Q and Q. Henceforth the first pair of stable cheapest finite element pair 
for the incompressible Stokes flows is given as 


<^1,0^ X ^0 with dimension 2iV* + Nq-2. (5) 


A second pair of stable cheapest finite element pair is obtained by enriching 
the velocity space by a globally one-dimensional Assnme that N is an even 
integer. Denote by the macro mesh snch that each macro rectangle 
consists of 2 X 2 rectangles Qjk,QjX+i^Q 3 +i,k^Q]+i,k+i, with (J,iX) = (j, fc). 
from 3"^ with J^K = 1,3, • • • , — 1. For each macro-element Q^x ^ 

define G DSSYg such that 

snpp{ipQM^) C qY, 

and the integral averages over the edges in vanish except 


/ 

dQj,k^dQj+i^k 




M ds = u, 

JK 


® ^pQM ds = —V. 

3 dQj^k+i<3dQjj^i^k+i 


where u denotes the unit outward normal vector to Qj^ on the edge dQji fl 
dQj+ie, i = k,k + 1. Introduce the following vector space of macro bubble 


functions: 


Span 






which is a one-dimensional subspace 
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Table 1: Number of degrees of freedom for different pairs (velocity/pressure) 


N 

^c,h o)dc,h 

<^ 2,0 ^ “*^ 1.0 

-^ 2.0 ^ "^0 

rotQi X 

*^ 1.0 ^ ^0 

2 ^ 

2178/289 

2178/255 

1088/255 

450/254 

2 ^ 

8450/1089 

8450/1023 

4224/1023 

192211022 

26 

33282/4225 

33282/4095 

16640/4095 

7938/4094 

2 " 

132098/16641 

132098/16383 

66048/16383 

32256/16382 


Table 2 : •^20 ^ *^20 ^ '^0 ’ *’otQi x and x stand for the two 

Taylor-Hood elements of type Q 2 x Qi and Q 2 x Pi, the Rannacher-Turek nonconforming 
quadrilateral rotated Qi x Pq element, and the nonconforming quadrilateral Pi x Pq element 
[ 28 ], respectively. 


of DSSYq. Then in enriched by addingdenoted by 

and hence, the second stable Stokes finite element pair is defined as follows: 

X ^g with dimension 2iV* + Nq. (6) 

The stability and optimal convergence properties of the two pairs of Stokes 
elements (5) and ( 6 ) are shown for the stationary Stokes equations in [28]. 

Comparing several other stable quadrilateral finite element pairs satisfying 
the inf-sup condition [13, 25, 31], the nonconforming element pairs (5) and ( 6 ) 
have the lowest degrees of freedom. Table 2 illustrates the degrees of freedom 
for different pairs, whose notations will be used throughout the paper. 

It is shown that both nonconforming finite element spaces (5) and ( 6 ) give 
exactly identical solutions for velocity fields but slight different pressure solu¬ 
tions whose differences in L^(f 2 )-norm are of order and thus both velocity 

and pressure are approximated with optimal convergence for Stokes flows. See 
[28, §4] for details. Due to this observation, we concentrate on the finite element 
pair (5) for approximating the cavity flow. 

3.4- Treatment of nonhomogeneous boundary condition 

In order to deal with the Dirichlet boundary values of cavity flows, the 
open boundary part (top boundary) is modified with two additional DSSY- 
type elements located at the two top corners, Qin and Qnn- Notice that for 
j = 1, • • ■ ,iV —I, has value I at the midpoints — I) and {{j+^)h, 1) 
and 0 at the other midpoints on the top boundary, one sees that 



assigns the vector value (1,0)* at the midpoints {{l + ^)h, I), • • • , {{N—l — ^)h, 1) 
and ( 5 , 0 )* at the midpoints (^h, 1) and (I — ^h,l), respectively. Denote the 
DSSY basis functions whose supports are the top corner elements as follows: 

= (Sf'' o fa-g and = (gf'' o . 
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both of which have values 1 at the top midpoints 1 ) and (1 — ^h, 1 ), re¬ 
spectively, and 0 at the other midpoints of the two elements. Summarizing the 
above, the approximate nonconforming finite element solution with the Dirich- 
let boundary data for the lid-driven cavity flow is approximated by Uh of the 
form 

Uh = uo,/i -b Ub^h, where (7) 

AT-l 

4>3,N + '4’2,TL + 1p2,TR ■ 

i=l 

We are now in a position to define a discrete weak formulation of (1) to find 
(uo,?i,p/i) G X ^1} such that 

) “b 5 ^ 0 ,/i 5 ? ^b,hi ^h) ~b Cfi ( 1 I 5 /^, , V/^) 

+bh(vh,Ph) = -ah(ub,h,Vh) - Ch(ub,h;ui,,h,Vh) Vv/^ G (8a) 

bh(uo,h,qh) = -bh(ub,h,gh) Vqh e ^ 0 , ( 8 b) 

where 

ah(uh,Vh) = ly / Vuh'.Vvhdx, bh(vh,qh) = - / (V ■ Vh)qh dx, 

Ch(wh;uh,Vh) = / (w/, • V)u/, • v/i dx. 

The nonlinear equations ( 8 ) can be approximated by the Picard iteration method 
[14, 16, 27]. With an initial guess (uQ°l,Pb^^) G ^" 0 ^ x ^q, define the Picard 

iterates (ug^l,p^f^^) G x for fc = 1, 2, • • • , solving the following Oseen 

problem: 

v^) -b bb(vh,Pb^) 

= -ah(ub^h,Vh) - + Ub,h;Ub,h,Vh) Vvh G (9a) 

bh(ulj^l,qh) =-bh(ub,h,qh) Vqh€^^. (9b) 

The Picard iterates (Ug^l,pj^^^) k>i are shown to converge at a linear order to 
the solution (uo,h,Ph) of ( 8 ) in [27]. One may of course use the Newton iterates 
which converge quadratically with sufficiently close initial guesses to the exact 
solution as described in [14, 16, 27]. 

4. Accuracy of solutions 

Previous studies validated their numerical solutions by comparing their nu¬ 
merical results with benchmark solutions in the literature, for example, [ 6 ] and 


Uo,/i 


Af-l 

= E 

i,fc=i 


Vjk 


^jk-i ^b,h — 




[19]. According to Erturk et al. [18], there are many different numerical proce¬ 
dures for the lid-driven cavity flow problem which yield very similar numerical 
results in the case of Re < 1000, however, their numerical solutions start to 
deviate from each other as the Reynolds number increases. 

Hence, in order to claim some sort of superiority of our nonconforming 
method over the other existing methods, we will not only compare our numer¬ 
ical results with those in the literature, but also show some other assessments 
for the accuracy of the numerical solution. 


4-1- Volumetric flow rate 

Aydin and Fenner [5] suggested a measurement of the accuracy of numerical 
solutions. They computed the net volumetric flow rate, Q, passing through a 
vertical line and a horizontal line to check the continuity of the fluid. Denote 
u = {u,v), and let Qu,c and Qv,c be the volumetric flow rate passing through 
a vertical line x = c and a horizontal line y = c, respectively. The volumetric 
flow rate values, Qu,c and Qv^c can be computed by 


Qu^C - 


u{c,y) dy 


Qv^c — 


v{x, c) dx 


( 10 ) 


4-2. Compatibility condition for the stream function tp 

We can also use the compatibility condition for the stream function ip for 
the assessment of the accuracy of numerical solutions. Using the expressions of 
the vorticity ui as the two-dimensional curl of the velocity: a; = V x u, and the 
velocity field as the two-dimensional curl of the stream function u = V x '0, one 
has the Neumann boundary value problem for ip as follows: 


— A'0 = Lo in D, 


with 

1 , 
1 , 
1 , 
1 . 

A compatibility condition, combined with (2), yields 


dtp 

dn 


= 


V for a; = 0, 

—V for a; = 1, 


/ oj dx = — 

Ja 



dtp 

dn 



( 11 ) 


( 12 ) 


One can compute oj dx by using the numerical solution u/i, and compare to 
check the accuracy of the numerical approximation. 
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4-3. Incompressibility condition 

Since the pointwise incompressible condition V- u = 0 should hold pointwise, 
the smallness of 


max / V ■ dx 


(13) 


is a good indicator to check numerical accuracy. This implies that (13) of the 
numerical solution should be close to zero. 

Invoking (7), and observing that 



one sees the following simplihcation: 

» N-l 

Y. / V • u„ dx = ^ 




'Q 


j,k=l • 


0(1)jk 0(1)jk 


N-l 






N 


i=i 


Ox 


dx. (15) 


Recall that (ftjk is piecewise linear, and hence its derivative is constant on each 
Qim- indeed, 




V(j)jk = < 




(-1,-1)* 


on Qjk, 
on Qj-\-i^kj 
on 

on Qj^k-\-i 1 


(16) 


for j = 1, • • • , iV - 1, fc = 1, • • • , iV. Set g/i = Effc=i CjkXQju ^ for a general 
piecewise constant element. By exploiting \Qf.m\ = b?, from (8b), (14), and (16) 
it then follows that 


N-l ^ N-l 

/ (V • \ih)qh dx = h ^ ) [^jkiCjk — Ci+i,fc+i) + Vjk{Cj+i,k ~ Ci,fe+i)] + ^ ( CjN 

j,k=i i=i 


As a basis for choose the union of the basis functions of q and q. 
For each j, fc = 1, 3, 5, • • • , A^ - 1, with qn = XQjk - XQn.n S and qt = 

XQj+i,k+i - XQn,n e foo™ (^b), (14) that 


V • dx = 


V • u/i dx = 


V • u?i dx. 


Qjk ^ Q ^ QNN 

Similarly, for j,k = 1,3,5, ■■■, N - 1, with qn = XQj+i,k - XQn-i,n € ^b,o 
Q h = XQj,k+i - XQn-i.n G ^b,o one concludes from (8b), (14) that 


(17) 

and 



u/i dx 



Uh dx 


V • dx. 


' Qn- 


1,N 


(18) 


10 








Setting 7 i = V • u/j dx and 72 = 1 n ^ obtains from 

(17) and (18) that 



However, using the Divergence Theorem piecewise for each Q G we have 


N 


/ V • u?i dx = y] / U-Uh ds = 'Y^ 


• uji ds = 0 . 






JdQjNn{y=l} 

Hence, ^(71 + 72 ) = 0 and therefore 

/ V • u,j dx = 71 VQjfc e 
JQjk 

( V • u,j dx = -71 MQjk e ^h- 

Qjk 


( 20 a) 

( 20 b) 


In order to compute 71 exactly, we sum ■ uji dx over all the red-type 

rectangles Qjk G invoking the form of u;i given in (7). Observing that for 
each interior vertex Vjk there are two rectangles in which share only the 
vertex: these two rectangles can be either the pair {Qjk,Qj+i,k+i) or the pair 
iQj+i,k,Qj,k+i)- Then the integrals / V • dx over those pairs cancel 

each other due to the Divergence Theorem or direct integrations. Recalling that 
the DSSY-type basis function parts on the two top corners do not contribute 
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anything for the integration of divergence, we see that 


/ V-u/,dx = y] / V • (uo^/i + dx 






Qjk 


= E 


V ■ U {,.,1 dx 


Qjk 


r /i\ 

E / (n) 


QjN(iSrR''Qi^ 


N/2 


i=i 


= E 


j = l JQ2j,N 
N/2 


V • ( Q ) {(t^ 2 j-l,N + 4>2j,N) dx 


= E 


— I J dQ2j,N 




{4’2j-l,N + 4>2j,N) ds 


N—1 / -i \ j 7 

'■eR) 


i=i 


( 21 ) 


A combination of (20) and (21) shows that 


V • Uft, dx 


Qjk 


= \/Qjk e ^h- 


( 22 ) 


We summarize the above results as in the following theorem: 

Theorem 4.1. Let ujj be in the form (7) and (uQ^h,Ph) Q x /^q fulfills 

(8b). Then (22) holds. Moreover, the signature of the integral over Qjk’s are 
alternating. 


5. Numerical results 

We have computed the steady state solutions of lid driven cavity flow from 
Re = 100 to Re = 5000 by using the Picard iteration method with the termina¬ 
tion condition: 


f - vAvl^^I - Nvl^^I - 

g — 

where A, B, and N denotes the matrices for the discrete Laplacian, divergence, 
and convection, respectively. First, notice that our proposed nonconforming 
finite element method for finding {uQji,Ph) G x fulfilling (8) does not 

modify the discontinuities at the top corners. However, the usual conforming 
hnite element methods require suitable modifications. For instance, for the 
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Table 3: Values used to plot the contours of the stream function and the vorticity 


Contours 

Values 





Stream function 

-0.1175, - 

0.1150, 

-0.11, -0.1, -0.09, 

-0.07, -0.05, ■ 

-0.03, -0.01, 


-l.OE-04, 

-l.OE-05, -l.OE-07, -l.OE- 

-10, l.OE-08, 

l.OE-07, 


l.OE-06, 

l.OE-05, 

5.0E-05, l.OE-04 

, 2.5E-04, 5.0E-04, 


l.OE-03, 

1.5E-03, 

3.0E-03 



Vorticity 

-5.0, -4.0, 

, -3.0, -2 

.0, -1.0, -0.5, 0.0, 

0.5, 1.0, 2.0, 

3.0, 4.0, 5.0 


■^2,0 X element method, the two popular cavity boundary conditions are 

used: the watertight cavity boundary condition 


g = 


( 1 , 0 )*, if 0 < a; < 1 and y = 1 , 
0 , elsewhere on d^l, 


(24) 


and the leaky cavity boundary condition 


g = 


( 1 , 0 )*, if 0 < a; < 1 and y = 1 , 
0 , elsewhere on d^l, 


(25) 


respectively. Notice that both conditions (24) and (25) satisfy (12). Our own 
FORTRAN and MATLAB codes were developed to implement the ^"g^ x 

element method while the IFISS S/W[l] was used to implement the ^ 2,0 ^ =^ 1 , 0 ^ 
element method. 

For the case of Re = 1000, we present in Fig. 5 the u-velocity profiles along 
the line x = 0.5 and the u-velocity profiles along the line y = 0.5 computed 
by using the x element with the boundary condition (3) and com¬ 

pare our results with those by Botella and Peyret [ 8 ] , by Bruneau and Saad [9] , 
and by Guermond and Minev [24]. In each case, our velocity profiles show a 
good agreement with the reference solutions. Recall that the solutions obtained 
Botella and Peyret used a spectral method on 160 x 160 spectral nodes, and 
Bruneau and Saad used a second-order finite difference schemes on the uni¬ 
form 1024 X 1024 nodes, while Guermond and Minev used a massively parallel 
computation combining their new direction splitting algorithm and the MAC 
central finite difference scheme on 5000 x 5000 nodes. 

For the range of Re = 100,400,1000,2500,3200, and 5000, Figs. 2 and 3 
show the w-velocity profiles along the line a; = 0.5 and the a;-velocity profiles 
along the line y = 0.5 computed by using the .^”g^ x J^g element with the 
boundary condition (3) and comparison results with those by Erturk et al. [18] 
and Ghia et al. [19]. In each case, our velocity profiles show a good agreement 
with their results. 

Although the velocity profiles in Fig. 5 and Figs. 2 and 3 seem to match 
quite well for Re = 1000, some of the actual numerical values differ in digits 
compared to those reported in [ 8 , 9, 24]. Hence, we compare the numerical 
values of the horizontal and vertical components of the velocity in Tables 4 
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Figure 1: Comparison of the vertical components of the rt-velocity along the segment x G [0,1], 
y = 1/2 and the horizontal components of the n-velocity along the segment x G [0,1], j/ = 1/2 
with Re = 1000 

Re=1000, NCj X Pj, (256x256) Re=1000, NCj x P", (256x256) 




an(i 5 with the reference solutions from [8, 9, 24]. Our numerical values, which 
were computed with 256 x 256 meshes with the lowest possible finite element 
X , match with the reference values mostly up to two digits, or with 
less than 1% errors; the numerical solutions, computed 512 x 512 meshes match 
with the reference values mostly up to three digits, or with less than 0.1% errors. 


Table 4: Comparison of the horizontal components of the velocity along the segment y G [0,1], 
X = 1/2 at Re = 1000 


y 

[8] 

[9] 

[24] 

X 0^1/(256 X 256) 

X l^^(512 X 512) 

0.0000 

0.0000000 

0.00000 

0.0000000 

0.0000000 

0.0000000 

0.0312 

-0.2279225 

NA 

-0.2279177 

-0.2274204 

-0.2276650 

0.0391 

-0.2936869 

-0.29330 

-0.2936814 

-0.2930076 

-0.2933552 

0.0469 

-0.3553213 

NA 

-0.3553154 

-0.3545665 

-0.3549485 

0.0547 

-0.4103754 

-0.41018 

-0.4103691 

-0.4096654 

-0.4100002 

0.0937 

-0.5264392 

NA 

-0.5264320 

-0.5271749 

-0.5264518 

0.1406 

-0.4264545 

-0.42645 

-0.4264492 

-0.4276315 

-0.4265356 

0.1953 

-0.3202137 

NA 

-0.3202068 

-0.3209943 

-0.3200577 

0.5000 

0.0257995 

0.02580 

0.0257987 

0.0256839 

0.0257175 

0.7656 

0.3253592 

NA 

0.3253529 

0.3259697 

0.3252217 

0.7734 

0.3339924 

0.33398 

0.3339860 

0.3346373 

0.3338694 

0.8437 

0.3769189 

NA 

0.3769119 

0.3778450 

0.3769140 

0.9062 

0.3330442 

0.33290 

0.3330381 

0.3339829 

0.3331021 

0.9219 

0.3099097 

NA 

0.3099041 

0.3108006 

0.3099725 

0.9297 

0.2962703 

0.29622 

0.2962650 

0.2971221 

0.2963312 

0.9375 

0.2807056 

NA 

0.2807005 

0.2815029 

0.2807605 

1.0000 

0.0000000 

0.00000 

0.0000000 

0.0000000 

0.0000000 


The computed streamlines are presented in Fig. 4. One can observe count- 
rotating secondary vortices at the bottom left and right corners of the square 
cavity. Bottom left and right vortices grow in size as Reynolds number increases 
and the secondary vortex at the top left corner of the square cavity develops as 
Reynolds number increases. 

The vorticity contours are presented in Fig. 5. We observe that the gradient 


14 














Figure 2: Profiles of u-velocity along the line 
with unregularized boundary condition 


Re=100, NCq (256x256) 



Re=1000, NC^ xP^i (256x256) 



Re=3200, NCj xpj, (256x256) 



= 0.5 computed by using the stable x 

Re=400, NC[| xP^i (256x256) 



Re=2500, NCq xP^^ (256x256) 



Re=5000, NCq xP^ (256x256) 
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Figure 3: Profiles of u-velocity along the line y = 0.5 computed by using the x 

with unregularized boundary condition 

Re=100, NCq (256x256) Re=400, NC[| xP^, (256x256) 
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Table 5: Comparison of the vertical components of the velocity along the segment x S [0,1], 
y = 1/2 at Re = 1000 


X 

[8] 

[9] 

[24] 

^" 0 ’'* X ^p{256 X 256) 

X g^l;{512 X 512) 

1.0000 

-1.0000000 

-1.00000 

-1.0000000 

-1.0000000 

-1.0000000 

0.9766 

-0.6644227 

NA 

-0.6644194 

-0.6666343 

-0.6648562 

0.9688 

-0.5808359 

-0.58031 

-0.5808318 

-0.5831751 

-0.5812660 

0.9609 

-0.5169277 

NA 

-0.5169214 

-0.5190905 

-0.5172781 

0.9531 

-0.4723329 

-0.47239 

-0.4723260 

-0.4741970 

-0.4725743 

0.8516 

-0.3372212 

NA 

-0.3372128 

-0.3380993 

-0.3370508 

0.7344 

-0.1886747 

-0.18861 

-0.1886680 

-0.1890994 

-0.1884232 

0.6172 

-0.0570178 

NA 

-0.0570151 

-0.0570951 

-0.0569011 

0.5000 

0.0620561 

0.06205 

0.0620535 

0.0622962 

-0.0619466 

0.4531 

0.1081999 

NA 

0.1081955 

0.1085611 

0.1080176 

0.2813 

0.2803696 

0.28040 

0.2803632 

0.2811184 

0.2802013 

0.1719 

0.3885691 

NA 

0.3885624 

0.3894565 

0.3885914 

0.1016 

0.3004561 

0.30029 

0.3004504 

0.3006758 

0.3004357 

0.0703 

0.2228955 

NA 

0.2228928 

0.2228075 

0.2228534 

0.0625 

0.2023300 

0.20227 

0.2023277 

0.2021815 

0.2022834 

0.0547 

0.1812881 

NA 

0.1812863 

0.1810885 

0.1812376 

0.0000 

0.0000000 

0.00000 

0.0000000 

0.0000000 

0.0000000 


in vorticity is negligible in the center of cavity and the region of very low gradient 
in vorticity grows as Reynolds number increases. 

In Table 6, we present the location of the center of the primary vortex, the 
stream function ijj, and vorticity oj at vortex center. These data are calculated 
for 100 < Re < 5000; for comparison, available data from the literatures are 
also given. The values of the stream function ip and vorticity ui are recorded 
at the center of meshes. The locations of primary vortices computed by using 
the .^"q^ X element differ from the other results by about 0.002 which is 
half the mesh size 1/h Ki 0.0039, Our numerical solutions computed by using 

both X and <^20 ^ ^ 10 ^ elements exhibit a good agreement with 

the literature data except in the case of the ^ 2 ’!o ^ ^ 10 ^ element with leaky 

cavity boundary condition (25) applied, ip and w for the x element 

with (3) are similar to those in the literature [18-20, 32]. In addition, Table 7 
summarizes data on the strengths and the locations of secondary vortices in 
the bottom left and right corners, and in the top left corner. We observe that 
secondary vortices appear stronger as the Reynolds number increases. 

In §4, we introduced the indicators for the accuracy of the numerical solu¬ 
tion. First, the volumetric flow rate values Qu,xc s-nd Qv,yc defined by (10) are 

shown in Table 8. for the <^” 0 ^ x element The values of Qu,xa and Qv,yc 
for the ^20 ^ element are much larger than those values at Qu,xc-h/ 2 : 

Qu,xc+h/ 2 , Qv,yc-h/ 2 , and Qy^y^+h /2 for the .^"g^ x ^1} element. Erturk et al. 
[18] calculated Q«,a;„and Qv,yc by using their solutions. The smallest values of 
Qu,xc = 4.5E-8 and Qv,yc = 1-34E-7 in [18] are lager than the largest values of 

Qu and Qy for the .^”g^ x .^g element. 

Table 9 shows the values of (12) and (13) for the x J^g and the <^20 ^ 

^j°g^ elements. Concerning the compatibility condition (12), the ,^"g^ x 
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Figure 5: Contours of vorticity computed by using the x element. 

Re=100, NCq xP^ (256x256) Re=400, NC^ (256x256) 
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Table 6: Computed primary vortex; the values of stream function (-0), vorticity (cj), and 
location (x,y). 


Rc 

FEM 

Grid 

'^min 

UJ 

[x, y) 

BC 

100 

*^2,0 

*^2,0 ’^1,0 

[19] 

[20] 

[32] 

256 X 256 

128 X 128 

128 X 128 

129 X 129 
128 X 128 

257 X 257 

-0.103531 

-0.103519 

-0.102872 

-0.103423 

-0.103435 

-0.103471 

3.16206 

3.18101 

3.15485 

3.16646 

3.1655 

(0.6152,0.7363) 

(0.6172,0.7383) 

(0.6172,0.7383) 

(0.6172,0.7344) 

(0.6172,0.7344) 

(0.6189,0.7400) 

(3) 

(24) 

(25) 

[20] 

[32] 

400 

£anc,h 

'^l.Q ^ ^0 

•^2,0 ”^1,0 

*^2,0 

[19] 

[20] 

[32] 

256 X 256 

128 X 128 

128 X 128 

257 X 257 
128 X 128 
257 X 257 

-0.114071 

-0.113990 

-0.111900 

-0.113909 

-0.113909 

-0.113897 

2.29821 

2.29476 

2.26041 

2.29469 

2.2950 

(0.5527,0.6035) 

(0.5547,0.6055) 

(0.5547,0.6055) 

(0.5547,0.6055) 

(0.5547,0.6094) 

(0.5536,0.6075) 

(3) 

(24) 

(25) 

[20] 

[32] 

1000 

V c:^h 

^1,0 ^ ^0 
£nc,h y c^dc,h 
“^2,0 "^1,0 

*^2,0 

[18] 

[19] 

[20] 

[32] 

256 X 256 

128 X 128 

128 X 128 

601 X 601 

257 X 257 
128 X 128 
257 X 257 

-0.119186 

-0.118941 

-0.115376 

-0.118781 

-0.117929 

-0.119173 

-0.118800 

2.07216 

2.06779 

2.00941 

2.06553 

2.04968 

2.0664 

(0.5293,0.5645) 

(0.5313,0.5664) 

(0.5313,0.5664) 

(0.5300,0.5650) 

(0.5313,0.5625) 

(0.5313,0.5625) 

(0.5335,0.5639) 

(3) 

(24) 

(25) 

[20] 

[32] 

2500 

cgtnc,h _ 

*^2,0 

V 

•^2,0 

[18] 

256 X 256 

128 X 128 

128 X 128 

601 X 601 

-0.122151 

-0.121492 

-0.115717 

-0.121035 

1.98912 

1.97645 

1.88476 

1.96968 

(0.5215,0.5449) 

(0.5195,0.5430) 

(0.5195,0.5430) 

(0.5200,0.5433) 

(3) 

(24) 

(25) 

3200 

£anc,h 

'^l.Q ^ ^0 

<5>c,/i ca^c,h 

*^2,0 “^1,0 

*^2,0 

[19] 

[20] 

[32] 

256 X 256 

128 X 128 

128 X 128 

257 X 257 
128 X 128 
257 X 257 

-0.122713 

-0.121860 

-0.115310 

-0.120377 

-0.121768 

-0.121628 

1.97778 

1.96186 

1.85833 

1.98860 

1.9593 

(0.5176,0.5410) 

(0.5195,0.5391) 

(0.5195,0.5430) 

(0.5165,0.5469) 

(0.5165,0.5352) 

(0.5201,0.5376) 

(3) 

(24) 

(25) 

[20] 

[32] 

5000 

(nc,h y (^dc,h 
*^2,0 

V 

*^2,0 

[18] 

[19] 

[20] 

[32] 

256 X 256 

128 X 128 

128 X 128 

601 X 601 

257 X 257 
128 X 128 
257 X 257 

-0.123658 

-0.122368 

-0.114120 

-0.121289 

-0.118966 

-0.121218 

-0.122050 

1.96650 

1.94277 

1.81321 

1.92660 

1.86016 

1.9392 

(0.5137,0.5371) 

(0.5156,0.5352) 

(0.5156,0.5352) 

(0.5150,0.5350) 

(0.5117,0.5352) 

(0.5156,0.5352) 

(0.5134,0.5376) 

(3) 

(24) 

(25) 

[20] 

[32] 


Table 7: Computed secondary vortices; the values of stream function ('i/j) and location {x^y) 
for the X element. 


Bottom left Bottom Right Top left 


Re 

t^max 

{x, y) 

l^max 

{x, y) 

t^max 

{x, y) 

100 

1.7368E-06 

(0.0332,0.0332) 

1.2597E-05 

(0.9434,0.0605) 

- 

- 

400 

1.4100E-05 

(0.0488,0.0488) 

6.4495E-04 

(0.8848,0.1230) 

- 

- 

1000 

2.3223E-04 

(0.0840,0.0762) 

1.7319E-03 

(0.8652,0.1113) 

- 

- 

2500 

9.2779E-04 

(0.0840,0.1113) 

2.6661E-03 

(0.8340,0.0918) 

3.3918E-04 

(0.0410,0.8887) 

3200 

1.1104E-03 

(0.0801,0.1191) 

2.S323E-03 

(0.8223,0.0840) 

7.0750E-04 

(0.0527,0.8965) 

5000 

1.3660E-03 

(0.0723,0.1387) 

3.0641E-03 

(0.8027,0.0723) 

1.4566e-03 

(0.0645,0.9082) 
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element and the Taylor-Hood element with the leaky cavity boundary condition 
(25) give precise values, while the Taylor-Hood element with the watertight 
cavity boundary condition (24) generates about 0.3% errors. An investigation 

of (13) shows that the numerical results obtained by using the x 

element are more accurate than those by the Taylor-Hood element. Moreover, 
the absolute values (13) for the x element are independent of Reynolds 

number and element Qjk due to Theorem 4.1. With the grid size 256 x 256, the 

absolute values (13) for the x element is given by 


max 

Qjk 



Uh dx 


1 

2563 


5.9605E-8, 


(26) 


while such values for the Taylor-Hood element with watertight and leaky cavity 
boundary conditions are given in Table 9. It should be stressed that the values 
obtained by the x element are smaller by a factor of four than those 

obtained by the Taylor-Hood element. 

At least judged by the three accuracy indicators, (10), (12), and (13), the 

numerical solutions by using the x element without any modification 

at the top corners are more accurate than those by using the Taylor-Hood 
element with modified boundary conditions (24) and (25). 


6. Conclusions 

The X .!3^g element is applied to solve the lid driven cavity problem 

with least modification at the two top corner element to deal with the jump 
discontinuities there using the DSSY element (of CDY element). 

The numerical solutions using g^ x element are compared with bench 
mark solutions and the horizontal and vertical components of the velocity at 
the center are correct up to mostly two and three digits if the mesh sizes are 
256 X 256 and 512 x 512, respectively. 

Numerical solutions were compared with those the conforming ^20 ^ ■^ 1 %^ 
element (Taylor-Hood element) with leaky and watertight cavity boundary con¬ 
ditions. Three indicators for accuracy of the numerical solution have been com¬ 
pared. (1) The incompressibility condition (2) The compatibility condition (3) 
with the Neumann boundary condition are used to check the accuracy of the 
numerical solutions. 

Our numerical solutions satisfy the incompressibility and compatibility con¬ 
dition precisely. Numerical results computed by using the x element 

show the best results in terms of satisfying incompressibility and compatibility 
conditions, and volumetric flow rates. 
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Table 8: Volumetric flow rates along the vertical and horizontal lines through the geometric 
center of the cavity, {xc,yc), by using the x and ^ 2 ’g ^ element 


Re 

FEM Grid Qu,., Qv,y. BC 

100 

X 256 X 256 1.9039e-16 1.2514e-13 (3) 

<^2,0 X 128 X 128 9.3009E-06 6.5662E-08 (24) 

^2.0 X 128 X 128 1.3114E-03 9.7804E-08 (25) 

400 

X 256 X 256 2.1554e-16 1.3347e-13 (3) 

-^2,0 X 128 X 128 1.4876E-05 1.2495E-06 (24) 

•^2,0 X 128 X 128 1.3170E-03 1.1132E-06 (25) 

1000 

X 256 X 256 3.5996e-17 1.1037e-14 (3) 

<^2,0 X 128 X 128 2.4097E-05 2.8794E-06 (24) 

•^2,0 X 128 X 128 1.3264E-03 2.5407E-06 (25) 

2500 

X 256 X 256 2.4373e-16 1.5280e-13 (3) 

<^2,0 X 128 X 128 4.0694E-05 5.7553E-06 (24) 

X 128 X 128 1.3431E-03 4.8544E-06 (25) 

3200 

X 256 X 256 2.1814e-16 5.1092e-14 (3) 

^2,0 X 128 X 128 4.6986E-05 7.0223E-06 (24) 

X .SZo^ 128 X 128 1.3494E-03 5.8405E-06 (25) 

5000 

-^Zo^ X 256 X 256 3.5562e-16 1.2311e-13 (3) 

■^2,0 X .^Zo^ 128 X 128 6.1055E-05 1.0206E-05 (24) 

X "^ZZ^ 128 X 128 1.3634E-03 8.2691E-06 (25) 
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Table 9: Compatibility (12) and incompressibility conditions (13) for the x and 

•^2 0 ^ elements. 


Re 

FEM 

Grid 

IJoW dx + l | 

(13) 

BC 


<::anc,h 

*^ 1,0 ^ ^0 

256 X 256 

2.8866e-15 

5.9605E-08 

( 3 ) 

100 

^c,h ^dc,h 

0 ^ 2,0 ^ *^ 1,0 

128 X 128 

2.6042e-03 

6.1596E-04 

(24) 


0C,h ^dc,h 

•^ 2,0 ^ *^ 1,0 

128 X 128 

1.1102e-15 

3.3407E-04 

(25) 


*^10 ^ "^0 

256 X 256 

2.2204e-16 

5.9605E-08 

( 3 ) 

400 

^c,k o)dc,h 

0 ^ 2,0 ^ *^ 1,0 

128 X 128 

2.6042e-03 

6.6730E-04 

(24) 


^»c,h 0)dc,h 

*^ 2.0 ^ *^ 1,0 

128 X 128 

4.7740e-15 

3.9730E-04 

(25) 


*^ 1,0 ^ "^0 

256 X 256 

2.6645e-15 

5.9605E-08 

( 3 ) 

1000 

^c,/i (j)dc,h 

0 ^ 2,0 ^ *^ 1,0 

128 X 128 

2.6042e-03 

7.2274E-04 

(24) 


^c,h (^dc,h 

*^ 2,0 ^ *^ 1,0 

128 X 128 

1.5543e-15 

5.0746E-04 

(25) 


ca'>T'C,h 

*^10 ^ "^0 

256 X 256 

1.1102e-15 

5.9605E-08 

( 3 ) 

2500 

^c,/i ^:T>dc,^ 

0 ^ 2,0 ^ *^ 1,0 

128 X 128 

2.6042e-03 

1.1836E-03 

(24) 


^c,h (^dc,h 

*^ 2.0 ^ *^ 1,0 

128 X 128 

7.7716e-15 

5.9441E-04 

(25) 


^10 ^ ^0 

256 X 256 

3.9968e-15 

5.9605E-08 

( 3 ) 

3200 

^c,h f^dc,h 

0 ^ 2,0 ^ *^ 1,0 

128 X 128 

2.6042e-03 

1.3685E-03 

(24) 


^c,h a)dc.h 

*^ 2.0 ^ *^ 1,0 

128 X 128 

1.5543e-15 

6.0657E-04 

(25) 


^10 ^ ^0 

256 X 256 

1.4433e-15 

5.9605E-08 

( 3 ) 

5000 

^c,h o)dc,h 

0 ^ 2,0 ^ *^ 1,0 

128 X 128 

2.6042e-03 

1.6240E-03 

(24) 


^c,h ^dc,h 

*^ 2.0 ^ *^ 1,0 

128 X 128 

4.2188e-15 

6.7909E-04 

(25) 
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